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Causal viscous hydro-dynamic fits to experimental data for pion and kaon transverse mo- 
mentum spectra from central Au+Au collisions at ^smn — 200 GeV are presented. Starting 
the hydrodynamic evolution at 1 fm/c and using small values for the relaxation time, rea- 
sonable fits up to moderate ratios rj/s ~ 0.4 can be obtained. It is found that a percentage 
£ — ' of roughly 50i]/s to 75r//s of the final meson multiplicity is due to viscous entropy produc- 

tion. Finally, it is shown that with increasing viscosity, the ratio of HBT radii i?out/-Rsidc 
approaches and eventually matches the experimental data. 

x> 

V : I. INTRODUCTION 



Extracting a value for the shear viscosity r\ out of data from the ongoing heavy-ion collision 
program at the Relativistic Heavy-Ion Collider (RHIC) is a difficult but maybe rewarding issue. It 
is difficult because the simplest hydrodynamic theory that includes viscosity, Navier-Stokes theory, 
is known to have problems with causality and instabilities in the relativistic case In order to 
repair these problems, so-called second-order theories have been put forward by Israel, Stewart 
[2] and Liu, Miiller, Ruggeri 0]. Unfortunately, this formulation introduces at least one other a 
priori unconstrained parameter, a relaxation time, into the hydrodynamic framework. Also, the 
C**"- ■ resulting hydrodynamic equations are quite complicated, and not easily implemented using existing 

numerical schemes 1 . 

Extracting rj from experiment may be rewarding, however, as results for rj from weak-coupling 



=5 



r^i i — l 



calculations for QCD [ll|, [l_2| and strong-coupling calculations for N = 4 SYM [13|, [pj] differ by 



an order of magnitude. Thus the value extracted from experiment might offer a clue whether the 



quark-gluon plasma at RHIC is described better by weak or strong coupling techniques 2 , a question 
currently hotly debated in the physics community. 

On the practical level, determining ry (or more precisely the ratio of shear viscosity over en- 
tropy density, rj/s) from experiment translates to finding the value of r\js for which a viscous 
hydrodynamic model fits experimental data best. Unfortunately, there is a lot of freedom in any 
hydrodynamic model calculation, mostly because of the poorly constrained initial conditions, but 
also in the "final" or freeze-out conditions (see below). It turns out that results for central collisions 
alone are not sufficient to fix all the free parameters, and as a consequence only upper bounds on 
the ratio of shear viscosity over entropy can be given in this work. Moreover, it should be noted that 
besides shear viscosity also bulk viscosity and heat-conductivity will affect hydrodynamic model 
fits to experimental data 3 . Algorithms on how to include these have been suggested in la, 1 1 Qf | . 



In order to keep complexity to a minimum, in the following only the effect of shear viscosity is 
included. This work is organized as follows: in section HU the setup of causal viscous hydrody- 



1 In the context of heavy-ion collisions, solutions for causal viscous hydrodynamic theories without transverse or lon- 
gitudinal dynamics (only "Bjorken flow") have been obtained in Refs. while results including transverse 
flow can be found in Refs. [1, H, [hj. No results including elliptic flow exist up to date. 

2 An anomalous viscosity value from turbulent magnetic fields [l5[ due to plasma instabilities [l6| ] could blur this 
naive picture, however. 

3 Bulk viscosity, however, has recently been calculated in QCD using weak-coupling techniques and found to be 
negligible compared to shear viscosity [13]. 
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namics as well as the initial conditions and equation of state for heavy-ion collisions are briefly 
reviewed. In section 1111} the equations for calculating particle spectra and HBT radii in viscous 
hydrodynamics are given and results are compared to experimental data in section IIVI Section [V] 
contains a summary and the conclusions. 



II. SETUP 

Including only the effects of shear viscosity, the causal viscous hydrodynamic equations used in 
the following are given by 0] 

(e+p)D/ = - A^ValT" 7 + W v Du v , (1) 

De = -(e + p)V/ + \w v (V u u^) , (2) 

TnA£A£L>II Q/3 + IP" = ??( W) - 2r n n Q ^w ly 2 , (3) 

where e,p are the energy density and pressure, respectively. The flow four- velocity obeys 
u^u^ = 1 and IP 1 ' is the shear tensor that fulfills u l JJ^ 1 ' = = 11^ and characterizes the deviations 
due to viscosity in the energy momentum tensor, 

p = (£+ P )«Y-j)f + n^. (4) 

The remaining definitions are 

d„u v = d^u u + T v ail u a , D = u^, W = A^du, 
M» = g ^ - u » u ^ = A^A^i (dpu a - d a u p ) , 

(A^B v )=A^B v + A v B^-l^ v A a B a , (A^, B v ) = \ (A fl B u + A U B^) , (5) 

where rj£„ are the Christoffel symbols. The parameter rn is a relaxation time that in weakly 
coupled QCD can be related to rj and the pressure as 0, 0] rn — which translates to rn ^ -jf;. 
To test for the dependence of the results on rn, in the following also a somewhat smaller value 



rn - will be used 4 , which can be argued for independently 20]. Note that formally one 
recovers the relativistic Navier-Stokes equations from Eq. (1,2,3) in the limit rn — > 0. 

The algorithm to solve the above equations including several tests was outlined in detail in 
Ref. [l0( and is not repeated here for brevity. In this work, the equations are solved on a lattice 
with 512 sites and a lattice spacing of 0.25 GeV -1 . 



A. Initial Conditions and Equation of State 



The energy density at the hydrodynamic initialization time tq is assumed to be parameterized 
by the number density of wounded nucleons in a Glauber model [211 ]. 



e( r 0) r ) = const x nwN(^) 



(6) 



nwNin 



T A (r) 



2T A (r) 



dz 



1-1 



A 
Po 



i + eX p[(V^T^-i?o)/x]' 



4 Once a non-trivial equation of state is used, final results will differ depending to whether one defines rn — — 



or m 



In this work, the latter definition is adopted, but final results seem to differ only slightly when 



implementing the other choice. 
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where po is such that 2ir J rdr Ta(t) = A. For gold nuclei, A = 197, i?o = 6.4 fm, \ = 0.54 fm 
and for numerical reasons ^1 — rjT ^- r ^ j is replaced by an exponential. The nucleon-nucleon cross 
section at ^/snn = 200 GeV is assumed to be given by a = 40 mb. Other parameterization of the 
initial energy density (e.g. scaling by the number of binary collisions), in general result in stronger 
radial gradients 2l| and consequently faster buildup of transverse flow. Since viscosity in some 
sense mimics the presence of transverse flow, the least restrictive initial condition and therefore 
the most conservative bound on viscosity will come from the choice of Eq. ([6]). In the same spirit, 
the initial value of IT/ 41 ' is chosen to be zero everywhere (see also the discussion in Ref. [nil]). 

The constant in Eq. (|6|) is chosen such that the central energy density e(r = 0) corresponds to 
a predefined starting temperature To via the equation of state. Since lattice QCD seems to rule 
out a first or second order phase transition [22T ]. the semi-realistic equation of state of Laine and 
Schroder 23] is used in the following. This equation of state is calculated from a hadron resonance 
gas at low temperatures and high-order weak-coupling QCD result at high temperatures with a 
cross-over transition near T c ~ 175 MeV (for brevity, the reader is referred to [23J for details). 



III. PARTICLE SPECTRA AND HBT RADII IN VISCOUS HYDRODYNAMICS 

A. Particle Spectra 



In order to convert hydrodynamic quantities into experimentally interesting observables, the 
standard method of choice is the Cooper- Frye freeze-out prescription 24\ . Assuming isothermal 
freeze-out at the temperature Tf defines a freeze-out surface which is characterized J25J] by its 
normal vector cffi^, 



cosh (77), cos ( 



dr(Q 

dC 



, sin 



• ; .sinl,l//)) r (C)r(C)^dC##- (7) 



The surface is parameterized by (e [0, 1] such that r(0) = tq and t(1) corresponding to the time 
when the last fluid element has cooled down to the temperature Tj. The single particle spectra for 
a particle with four momentum j/ 4 = (E, p) and degeneracy d are then calculated as 



E 



d 3 N 



d 



d*p (2vr) s 



P^f 



11/' 



(8) 



where it is reminded that the viscous distribution function / is related to the ideal distribution 
fo(x) = [exp(x) ± 1] _1 and the shear tensor IP 1 ' as 0, 26] 



/ 



T 



fo 



T 



lT/o 



1 + 



p„p v W 



2T 2 {e+p) 



(9) 



where ± applies for fermions and bosons, respectively. For simplicity it is, however, convenient 
to approximate fo(x) ~ exp (— x) which does not seem to affect final results too much. In this 
case, one consequently also has to replace the expression [1 t- f(x)] in Eq. © by 1 and all but one 
integral in Eq. ([8]) can be evaluated analytically to give [l(| 



E- 



J 3 N 
d 3 p 



2d 
W) 2 Jo 



d( r(C) r(C) 



dr 
d( 



m_ L I (u r p ± /T)K 1 (u T m_ L /T) 



p±Ii(u r p ± /T)K (u T m ± /T) 



+ 



2pxjn\_vW T 
4T 2 (e + p) 



— m ± (K + K 2 )h - j-piKxih + I ) 
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+ 



4T 2 (e + p) 

4T 2 {e + p) 

v 2 p\W r 
4T 2 (e + p) 

4T 2 (e + p) 

m 2 ± m v 

~AT 2 {e+p) 



dv cLt 

dv 1 ClT 

— -m^SKx + K 3 )I - j^Pl(K + K 2 )I X 

dv cIt 1 

— m ± K 1 (I - I 2 ) - — p±K -(Ii - h) 

dr 



m ± K 1 (I 



dr 1 (v 



Ki)I -^p ± (K 2 



K )h 



(10) 



where I n (x) and K n (y) are modified Bessel functions that have arguments x = u r p±/T and y 
u T m±/T as denoted in the first two lines. 



B. HBT Radii 



Given two identical particles with four momenta p^ and p^, respectively, the coincidence prob- 
ability of measuring these two particles in a single event divided by the probability of the particles 
being uncorrelated defines the two-particle correlation function Cip^p^. Rewriting the corre- 
lation function as a function of the momentum difference = P\ — p 2 and average momentum 
= 2 {pi + p 2 ) and assuming chaoticity and large size of the emitting particle source, the cor- 
relation function may be written as 27] 



1 + 



E\E 2 



d 6 N 
a! 3 pi<i 3 P2 



TP dgjV jp d 3 N 



(11) 



P2 



where ^J^-, E 2 -j^j- are the single particle spectra defined earlier and 



E\E 2 



d 6 N 



d 3 pid 3 p 2 



(2vr) s 



(12) 



If furthermore boost-invariance and rotational symmetry around the longitudinal axis is imposed, 
one can choose the average transverse momentum as k = (k x ,k y ,k z ) = (fc, 0,0) and decompose 
q into so-called "out", "side" and "long" components, q ou t = (<?out, 0, 0), q s jd e = (0, </ S ide> 0), 
qiong = (0, 0, giong)) respectively. The correlation function C(p±,p 2 ) for these components then 
also can be used to define the three HBT-radii via the Bertsch-Pratt parameterization [281 . |29| 



C(q, k) ~ 1 + exp -R 2 out (k)ql ut - R 2 idc (k) Qsidc ~ RfonMql 



ong 



(13) 



Even though the correlation function C obtained through Eq. (11 If) will in general not have the above 
Gaussian form, for simplicity the HBT radii are determined as R = \/ln2/g* where C(q*, k) = 1.5, 
as proposed in Ref. [s| (see also Ref. [25]] ) . 



With the freeze-out surface taking the form 25] 

(S*, E 2 , E», £ z ) = (r(C) coshfa), r(() cos(0),r(C) sin(0), r(Q sinhfa)) (14) 

one can again do some of the integrals in Eq. (|12|) analytically. Specifically, it is found that for 
i?i ong , the result for / s dS^i^^exp [iT,^]/ (K^u^) becomes an expression similar to Eq. (fTUj) . 
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FIG. 1: Viscous hydrodynamic fits of pion spectra to experimental data from PHENIX (circles) and STAR 
(squares) [33l [3^|. for relaxation time values Tn = (left) and rn = (right). Also shown are results 
for kaons and protons, scaled by 0.1 and 0.01, respectively. (The hydrodynamic results and the STAR data 
include weak decays while the PHENIX data does not). See text for details. 



but with all Bessel K functions replaced by 



K n (x) ->■ K n 



cosh (n rj) cos [rgi on g sinh rj\ exp [—x cosh rj\ . 



For iZside; one nas t° replace all Bessel I functions in Eq. (jlOp by 



(15) 



In{x) — ► I n = / — cos (n (f>) cos [rg si( j e sin((^>)] exp [x cos 

JO T 

and finally, for R ou t, the arguments of the Bessel I and K functions have to be replaced by 



x -> x = — it (Ei - E 2 ) 



u m± 

y -^y = * r <iout- 



(16) 



(17) 



(Note that in this case the modulus in Eq. (|12p is important.) Parts of these simple relations have 
been found in 0, 25]. 



IV. RESULTS 



A. Meson Spectra 



Spectra of pions and kaons have contributions from two sources: firstly, there are direct contri- 
butions for these particles at freeze-out, which are calculated using Eq. (fTU|) . Secondly, there are 
contributions from unstable hadrons and hadron resonances that can decay into pions and kaons. 
The spectra of these unstable particles with masses up to 2 GeV are also calculated at freeze-out 
via Eq. (|10p and then all possible two- and three-body decays that contribute to the stable particles 
of interest are taken into account using the decay routine from the AZHYDRO [i3] code, based on 
Refs. HILEI]. 

The initial central temperature To and the freeze-out temperature Tf are adjusted such that 
both the normalization and slope of the resulting pion spectrum is in reasonable agreement with 
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dA'A'.visc / dN K idcal 
dy 1 dy 


rj/s = 0.08 


1.06 


1.06 


T)/s = 0.16, r n = f f 


1.12 


1.12 


V /s = 0.16, r n = ^ 


1.12 


1.12 


jj/a - 0.24, r n = f 6 


1.15 


1.15 


ry/s = 0.24 r n = § J# 


1.18 


1.19 


77/s = 0.32 


1.23 


1.23 


rj/s = 0.40 


1.28 


1.28 



TABLE I: Fraction of multiplicity due to viscous entropy production: shown are the ratios of total multi- 
plicities for pions and kaons calculated in viscous and ideal hydrodynamics, for identical initial/freeze-out 
conditions. 

experimental data J32). l34| . This is the procedure adopted in ideal hydrodynamics [3H], and - as 
anticipated in Ref. [y, ll0(] - can be carried through also for non-vanishing rj/s, as is shown in Fig.[TJ 
Usually, in addition the initialization time To is allowed to vary in order to obtain a best fit also 
of the elliptic flow experimental data. Since in this work only central collisions are studied, the 
choice to = 1 fm/c is adopted, but results should not depend strongly on this choice. 

In Fig. [H results for the spectra of pions, kaons and protons for various hydrodynamic runs 
with different values of shear viscosity are shown together with experimental data from PHENIX 
[33| and STAR [3] experiments for the most central 5% of Au+Au collisions at y / J/v r /v = 200 GeV. 
It should be noted that no chemical potential is included in the equation of state, and therefore a 
distinction between particles and anti-particles is not possible. As a consequence, the spectrum of 
protons cannot be expected to match the experimental data, but is included in order to demonstrate 
that it is at least close to the experimental result. 

In general, increasing the value ofrj/s and leaving the initial and freeze-out conditions unchanged 
tends to make resulting spectra flatter 0, 0, 03, [2a] and hence in some sense mimics increasingly 
stronger transverse flow. Starting from initial/freeze-out conditions for which ideal hydrodynamics 
fits the experimental spectra and smoothly "turning on" viscosity, one has to counteract the effect 
of rj/s by reducing the buildup of transverse flow to keep the spectra in agreement with the 
experimental data. Both decreasing the initial central temperature and increasing the freeze- 
out temperature reduces the amount of hydrodynamic transverse flow at freeze-out. The former 
significantly affects the total pion multiplicity whereas the latter does not, offering a convenient 
way of decreasing total transverse flow while keeping the overall meson multiplicity close to data. 
However, since the concept of the Cooper-Frye freeze-out mechanism including the effects from 
resonance decays probably does not make sense at temperatures that are too high (far above 
T c ~ 175 MeV), an upper limit TV < 200 MeV is imposed. This limiting TV is much larger than what 
is typically used in ideal hydrodynamic model fits. However, in ideal hydrodynamics the system 
interactions by definition always keep the system in perfect thermal equilibrium, while the presence 
of viscosity means that interactions are not as efficient, allowing for departures from equilibrium. As 
a consequence, an earlier freezing-out and thus a higher TV for viscous hydrodynamics as compared 
to ideal hydrodynamics is to be expected. 

In Fig. [TJ results are shown for a relaxation time value of rn = (left) and r n = ^ ^ (right) 
and various rj/s. For rn = ^ the values of To used to generate the figure are To = 0.34, 0.33, 0.33 
GeV for rj/s = 0.08,0.16,0.24. For r n = ffi the values are T = 0.34,0.32,0.31 GeV for rj/s = 
0.16, 0.32, 0.40. From this figure, it can be seen that for the larger value of rn, the pion (and kaon) 
spectra in viscous hydrodynamics can be made to agree reasonably well with experimental data 
up to rj/s ~ 0.25 by mainly changing the freeze-out temperature TV. At higher values of rj/s, 
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FIG. 2: Pion HBT Radii from hydrodynamics compared to data from the STAR experiment [361 ]. Shown 
are the ratios R ou t/R s ide (left) and -Riong/-Rsidc (right) for ideal and viscous hydrodynamics (with rn = -% 
except for results indicated by arrows where rn 
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the pion spectrum becomes either too flat (when leaving Tq,Tj unchanged), the pion multiplicity 
becomes too low (when To is lowered such that the slope matches the data) or one has to choose 
an unreasonably large Tf > 200 MeV. 

The same is true when adopting a smaller value of rn, but it is possible to have reasonable 
agreement with experimental data up to rj/s ~ 0.4 (see Fig. [1]), about twice the bound from the 
larger value of rn- A more detailed study of the effect of changing rn is currently in progress 20]. 

For non-vanishing viscosity, the hydrodynamic evolution is no longer isentropic. Some part of 
the final multiplicity is therefore due to entropy production; in order to give some quantitative 
answers how much entropy is created, it is instructive to consider the ratio of meson multiplicities 
at non- vanishing rj/s and ideal hydrodynamics for identical initial/freeze-out conditions. Using 
the initial/freeze-out conditions for which the viscous calculations match the experimental pion 
spectrum, one finds the results shown in Tabled It seems plausible to infer from Table U that final 
multiplicities are increased by a fraction of ~ 0.75 rj/s with respect to ideal hydrodynamics due to 
viscous entropy production. Roughly, this translates to a percentage from 50 rj/s to 75 rj/s of the 
final multiplicity being due to viscous effects. 



B. HBT Radii 



Using the definition of the HBT radii from the previous section and the initial/freeze-out con- 
ditions for which the viscous hydrodynamic calculation matches the experimental pion spectrum, 
one can proceed to compare the viscous hydrodynamic results for the HBT radii to experimental 
data. This is done in Fig. [2J where the ratios R on t/R s ide arid R\ on g/ Rside for various values of rj/s 
are compared to data from STAR 36|. The ratio Ri ong /Rout is of special interest since in ideal hy- 
drodynamics it has been notoriously hard to obtain results that are not far above the experimental 



data (which sometimes is known as the "HBT puzzle"). It has been argued in Refs. [1, 26, 37] 
that non-vanishing viscosity may decrease this ratio and thus bring it closer to the data. Fig. [2] 
represents the first result of the ratio R\ on g/Rout for non- vanishing viscosity and initial/freeze-out 
conditions for which the experimental pion spectrum is matched at the same time. 

As can be seen from this figure, the ratio i? out /i2 s ide approaches the experimental data as 
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FIG. 3: Pion HBT Radii from hydrodynamics compared to data from the STAR experiment j36|. Shown are 
-Rout (left) and i? S idc (right) for ideal and viscous hydrodynamics (with m = except for results indicated 
by arrows where rn = f^)- See text for details. 



viscosity is increased. Indeed, for all but the lowest p± values the ratio is consistent with the 
data (within error bars) for the highest value of viscosity where the pion spectrum still matches 
the experimental data. Remarkable as this may be, it is unlikely that the presence of viscosity 
completely solves the HBT puzzle, the reason being that even though the ratio -Rout /-Rside in the 
viscous hydrodynamic calculations moves close to data, the absolute values R ut, Rside tend to 
be below the experimental values, as can be seen from Fig. [3j Put differently, the agreement of 
R ou t /Rside with data is achieved mainly by lowering i? ut, while i? s ide is hardly affected and always 
stays much below the experimental values. 

At present, it is not excluded that a match of the viscous hydrodynamic -R s ide with experimental 
data can be achieved when e. g. w hen changing also the hydrodynamic initialization time To, or the 
initial energy density profile [35]. However, the mismatch of i? s ide with data could also indicate 
that other effects, such as the replacement of the sharp hypersurface freeze-out by a more realistic 
continuous emission of particles 38|,l39|] have to be taken into account to obtain a proper description 
of the experimental HBT radii. 

Also shown in Fig. [2] is the ratio R\ ng/ Rside, compared to experimental data [36]. In this case, 
the viscous hydrodynamic result can be made to intersect the experimental results, but the too 
steep slope seems to prevent a nice match to data. Again the absolute values tend to be below the 
experimental data, also perhaps indicating the need for a more realistic freeze-out treatment. 



V. CONCLUSIONS 



Using a simple numerical code to solve the causal viscous hydrodynamic equations for the case 
of central heavy-ion collisions, it has been shown that both pion total multiplicities and spectral 
slopes for Au+Au collisions at ^/snn = 200 GeV can be matched to experimental data for moderate 
viscosities 5 . The bound on rj/s up to which this matching is possible depends on the value of the 
relaxation time (which is a priori a free parameter of the causal viscous hydrodynamic framework) 



5 The C++ code including results is available from 



http://hep.itp.tuwien.ac.at/' [paulrom / 
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and on the hydro-dynamic initialization time tq. Setting tq = 1 fm/c, one obtains rj/s < 0.25 for a 
value of the relaxation time rn that corresponds to the weak-coupling QCD result, and rj/s ^ 0.4 
for a somewhat smaller value of rn- For all the viscous hydrodynamic calculations that result in a 
pion spectrum matching the experimental data, it was found that roughly 50 r//s to 75 rj/s percent 
of the final multiplicity is due to viscous entropy production. 

For hydrodynamic parameters that allow to match the experimental pion spectrum, it was 
shown that the ratio of HBT radii i? ut/-Rside approaches and eventually matches the experimental 
result as r\j s is increased. The absolute values of i? ut, -Rside; however, tend to be much below the 
experimental values, suggesting that a more elaborate freeze-out treatment (or other effects, see 
e.g. Ref. [40]) might be necessary to achieve agreement. 

The above bound on the ratio n/s ^ 0.4 is interesting since it lies roughly half-way in between 



the weak-coupling QCD result [111. Il2l] and the strong coupling AT = 4 SYM result [13j, [lj]. If this 



bound was saturated, this could suggest that the quark-gluon plasma created at the highest RHIC 
energies is neither a very weakly nor a very strongly coupled plasma, but rather a border-case 
between these two. While this may not be the most elegant of possibilities that nature could have 
chosen, a not-weakly, not-strongly coupled plasma which behaves as a non-ideal fluid might be 
the most realistic assessment of the current status of RHIC data (see Ref. 41. 142I 43] for related 
conclusions) . 

To decide, and maybe extract a value of the ratio r]/s (and at the same time rn) from RHIC 
data, it seems necessary to extend this work to the treatment of non-central collisions, and thus 
viscous hydrodynamic results for elliptic flow. 
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